Critical temperature for the two-dimensional attractive Hubbard Model 
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The critical temperature for the attractive Hubbard model on a square lattice is determined from 
the analysis of two independent quantities, the helicity modulus, p 3 , and the pairing correlation 
function, P s . These quantities have been calculated through Quantum Monte Carlo simulations for 
lattices up to 18 x 18, and for several densities, in the intermediate-coupling regime. Imposing the 
universal-jump condition for an accurately calculated p s , together with thorough finite-size scaling 
analyses (in the spirit of the phenomenological renormalization group) of P 3 , suggests that T c is 
considerably higher than hitherto assumed. 
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The attractive Hubbard model[l],l2 has been success- 
fully used to elucidate a number of important and funda- 
mental issues in both conventional and high-temperature 
(cuprate) superconductivity. The nature of the crossover 
between BCS superconductivity (at weak coupling, or 
small on-site attraction) and Bose-Einstein condensa- 
tion of tightly bound pairs (strong coupling) has been 
shown to be smooth 0, The appearance of pre- 
formed pairs within a certain range of parameters in the 
normal phase, especially below a characteristic temper- 
ature, has been related to pseudogap behavior of high- 
temperature superconductors 0, [J . Further, this model 
allows one to introduce disorder on the fermionic degrees 
of freedom 0, Q and investigate the behavior near the 
quantum critical point of the two-dimensional insulator- 
superconductor transition; this provides an alternative 
to the dirty-boson picture 9] to discuss the universal 
conductivity [ljj. The attractive Hubbard model with 
a periodic modulation of U has been used to interpret 
superconductivity in layered structures 

A basic concern has run through many of these cal- 
culations, in particular those based on Quantum Monte 
Carlo (QMC) simulations. In two dimensions, there is a 
consensus that the early QMC phase diagram [TJ, |l3| - in 
the space of critical temperature, T c , electronic density, 
(n), and magnitude of the on-site attraction, \U\ - is qual- 
itatively correct. However, some serious quantitative dis- 
crepancies have emerged over the years, pointing towards 
higher critical temperatures; see, e.g., the Bogoliubov- 
Hartree-Fock (BHF) approach of Ref. Our purpose 
here is to examine the dependence of T c with (n), for 
fixed U, by resorting to a much wider (namely, larger 
system sizes and several electronic densities) set of QMC 
data, together with alternative procedures to locate the 
critical temperature. Establishing an accurate value for 
this most fundamental property of the model is impor- 
tant, especially as the physics of variants of the attractive 
Hubbard Hamiltonian is explored, and comparisons are 
made to the original system. 



The model is defined by the Hamiltonian 
U = - l E ( c Lcj CT + H.c.)- M ^(n iT +n. 

(U>,<r i 
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where c. ff (c icr ) destroys (creates) an electron with spin a 
on site i of a square lattice, (i, j) denotes nearest neighbor 
sites, n- lrT = c ia c ia , \U\ is the strength of the attractive 
interaction and fj, is the chemical potential. From now 
on, all energies are expressed in units of the hopping am- 
plitude, t, and we also set fcg = 1. 

At half filling (which corresponds to the particle-hole 
symmetric point, fj, = 0), the degeneracy of charge- 
density wave (CDW) and singlet superconducting (SS) 
correlations leads to a three-component order parameter 
[lfij : the transition temperature is therefore suppressed 
to zero. Away from half-filling, CDW correlations are 
suppressed and a finite temperature Kosterlitz-Thouless 
(KT) transition into a SS phase takes place (l^.ll^|: 
this phase has only algebraically decaying correlations 
for < T < T c . Further, close to half filling an exact 
mapping onto the two-dimensional Heisenberg antifer- 
romagnetic model in a magnetic field leads to 0, 
T c ~ — 27rJ/ln|l — (n)|, so that T c rises sharply from 
zero as one dopes away from (n) = 1. 

We start by employing the analysis of Ref. to new 
data for the SS pairing correlation function, 



P s = (AtA + AAt), 



with 



At = 



E c iT c U- 



For < T < T c , one expects 
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Figure 1: (Color online) P 3 as a function of /3 = 1/T for 
(n) = 0.5 and different lattice sizes, L. 



where r = |i — j|, and rj(T) increases monotonically be- 
tween r)(0) = and r)(T c ) = 1/4 

The hnite-size scaling behavior of P s is therefore ob- 
tained upon integration of T(r) over a two-dimensional 
system of linear dimension L. One then has |13j 



P s = L 2 ~^f(L/0, L > 1, T -> T+, (5) 



with 



£ ~ exp 



.4 



(T - T c )i/2 



(6) 



in the thermodynamic limit, one recovers P s ~ £ 7 / 4 . 
For completeness, one should mention that since ?y — > 
as T — > 0, the system displays long-range order in the 
ground state, so that a 'spin-wave scaling' is expected to 
hold fil: 



(7) 



where Ao is the superconducting gap function at zero 
temperature, and C is a |£7|-dependent constant. 

Similarly to Ref. Il3l here we use the determinant 
QMC algorithm [2(J to calculate P s . Typically our data 
have been obtained after 500 warming-up steps followed 
by 50,000 sweeps through the lattice. The discretized 
imaginary-time interval [20j was set to At = 0.125, which 
is small enough for the results not to depend on this 
choice in any significant way. 

In Fig. n we show raw data for P s as a function 
of (1/T), for fixed density, (n) = 0.5 and U = -4, 
and for different lattice sizes. The crossover between 
temperature- and size-limited reg imes is described by 
finite-size scaling (FSS) theory |2j| and appears as a lev- 
elling off of P s below a certain temperature for each sys- 
tem size. Before a more quantitative scaling analysis, we 



Figure 2: (Color online) Rescaled P 3 as a function of w = 
Lexp[-A/(T-T c ) 1/2 \ for (n) = 0.5 and different lattice sizes, 
L. The values for A and T c are the ones determined in Ref. 
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can already see a suggestion that T c is around 1/6 from 
the raw P s data. In general, at temperatures for which 
correlations are short ranged, a structure factor like P s is 
independent of lattice size. As T is decreased, the point 
at which the structure factor begins growing with lat- 
tice size signals the temperature at which the correlation 
length £ is becoming large (comparable to the lattice size 
L), thus providing a crude estimate of T c . The subse- 
quent plateau at low temperatures occurs when £ >> L. 
This crossover is contemplated by the FSS form, Eq. 
(5), which can be invoked to determine T c by plotting 
L- J / 4 P S as a function of of w = Lexp[-A/(T - T c ) 1/2 ], 
at a given U, for different system sizes, with T c and A 
being adjusted to give the best possible data collapse, 
as done in Ref. Figure [3 shows the resulting scaling 
plot, in which the values A = 0.4 and T c — 0.045 were 
determined in Ref. With our substantially increased 
amount of data points it becomes clear that the data col- 
lapse onto a single curve with the parameters A and T c of 
Ref. becomes rather unsatisfactory. We furthermore 
note that Eq. © is expected to hold only for t < 10~ 2 
(t = (T — T c )/T c ), see, e.g., Ref.|13 For the value of T c 
used in Fig. |3 only few data obtained for L = 4, 6, 8, 10 
m Ref. [HI satisfy this criterion. 

We therefore obtain new values of A and T c from our 
expanded data set. We disregard the data from the small- 
est system sizes, L = 4 (which, additionally, has a spe- 
cial topology, being equivalent toa2x2x2x2 four- 
dimensional lattice), L = 6, L = 8 and also L = 10; 
as we will see below, the SS pairing correlation function 
presents large finite size effects for these lattice sizes. Fur- 
thermore, we only include data points for temperatures 
T for which Eq. © is expected to hold (see above) . Fig. 
[3]clearly shows that, for the larger latices, our newly de- 
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Figure 3: (Color online) Same as Fig. [5] but with A and T c 

determined from the present data. Inset shows same system Figure 4: (Color online) Helicity modulus as a function of 
sizes, with A and T c from Ref. 0. temperature for (n) = 0.5 and different lattice sizes, L. The 

straight line corresponds to 2T/n. 



tcrmincd parameters A = 0.1 and critical temperature 
T c = 0.13 render a much better data collapse than the 
old parameters. 

The present analysis shows that the estimates of T c ob- 
tained in Ref. ^3 can be quantitatively quite unreliable. 
In our opinion, this is due to the fact that the finite-size 
behaviour of P s , Eqs. 101-®! follows from an analysis 
which is valid only for large enough lattice sizes, since it 
involves the binding-unbinding of rather large structures 
(vortices) in the KT transition 16]. Moreover, the pa- 
rameters A and T c that have to be found via data fitting 
both reside in an exponent, resulting in large uncertain- 
ties for the individual fitted parameters. Although the 
lattice sizes used in the present study may not be large 
enough to determine T c with high accuracy, our result is 
bound to be an improvement and in any case indicates 
that the actual T c may be much larger (by even a factor 
of 3) than believed sofar. 

This tendency towards higher critical temperatures ap- 
pears as well in a completely independent analysis, based 
on the behavior of the helicity modulus (HM). The latter 
is a measure of the response of the system in the ordered 
phase to a 'twist' of the order parameter j2^|, and can 
be expressed in terms of the current-current correlation 
functions as follows I2JI 
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[A L - A T ] , 



where D s is the superfluid weight, and 

A L = lim A xx (q Xl q y = 0,uj n = 0), 

and 

A T = lim A xx (q x = 0,q y ,U! n = 0), 

q v ^0 



(8) 



(9) 



(10) 



are, respectively, the limiting longitudinal and transverse 
responses, with 

A x Mco n ) = / dr e^e^ T A xx (£,T), (11) 



where uo n — 2nnT; 



A xx (e,T) = (j x (i,T) Jx (o,o)), 



(12) 



where 



Ulr) = e^EtW^ " 4 W ]e-^ (13) 

a 

is the x-component of the current density operator; see 
Ref. El for details. 

At the KT transition, the following universal-jump re- 
lation involving the helicity modulus holds |2oj|: 



7T 

Tc = - Ps 



(14) 



where p~ is the value of the helicity modulus just be- 
low the critical temperature. Thus we can obtain T c 
by plotting p s (T), and looking for the intercept with 
2T/t:. This procedure has been used before, with p s 
calculated within a Bogoliubov-Hartree-Fock (BHF) ap- 
proximation 0,H3; since transverse current-current cor- 
relations were neglected, p s is likely to have been overes- 
timated, and the ensuing T c 's may have been too high. 
Here we calculate both A L and A T by QMC simulations 
to obtain p s through Eq. (JSJ); a typical example of p s (T), 
for (n) — 0.5, is shown in Fig. 0] We see that finite- 
size effects are not too drastic, since all curves cross the 
straight line within a small range of temperatures; that 
is, from Fig.^Jwe can estimate T c = 0.14 ± 0.02. 
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Figure 5: (Color online) Logarithm-linear plot for the rescaled 
P s as a function of (3, for (n) = 0.5 and for different lattice 
sizes, L, symbols are the same as in Fig. 2. The inset shows a 
blow-up of the region centered about (3 = 7. No parameters 
are adjusted. 



In order to check the robustness of this method, we can 
extract T c from P s through a 'phcnomenological renor- 
malization group' (PRG) |27ll28j analysis, provided some 
subtleties peculiar to the KT transition are kept in mind. 
Since £ — > oo for all T < T c , Eq. JSJ implies that curves 
for L~ J / 4 P S (L, /?), when plotted as functions of f3, and 
for different L, should all merge for f3 > f3 c . Figure El 
shows that this characteristic feature only sets in for the 
largest system sizes, namely, L > 12, from which we can 
infer f} c — 7.5 ± 0.25; these error bars are somewhat ar- 
bitrary, and result from visual inspection. It should be 
stressed that this estimate for (3 C agrees remarkably well 
with the one obtained from the helicity modulus, indi- 
cating the robustness of both procedures to extract T c . 
Interestingly, we should notice that the curves for L = 6 
and 8 cross each other (as in an ordinary second order 
transition) at [3 ~ 7, which is very close to (3 C estimated 
from the larger systems. Therefore, within the context 
of PRG, for the smallest sizes a KT transition appears as 
an ordinary transition, only crossing over to the merging 
feature for the largest sizes. 

The critical temperature has been estimated for other 
electronic densities, (n) = 0.1 (HM and PRG), 0.3 
(PRG), 0.7 (PRG), and 0.875 (HM and PRG); all PRG 
plots display the crossing and merging tendency observed 
for (n) — 0.5. The resulting phase diagram is shown in 
Fig. El for comparison, we also plot the early QMC results 
pL3| , the parquet data from Luo and Bickers [2!j , and the 
estimates from the BHF approximation [l4| . While close 
to half filling all results (but BHF) are in fair agreement, 
for larger dopings agreement is only found between the 
results from PRG and those from the helicity modulus. 
The inescapable conclusion is that the critical temper- 



Figure 6: (Color online) Critical temperature as a function of 
band-filling, obtained by different methods. All lines through 
data points are guides to the eye. 

ature for the superconducting transition in the attrac- 
tive Hubbard model is actually higher than previously 
assumed. 

In summary, we have established very reliable esti- 
mates for the critical temperature of the square-lattice 
attractive Hubbard model. This is done by finding quan- 
titative agreement between entirely different procedures 
to extract T c from two independent correlation functions 
(both computed by determinant QMC). As a result, the 
critical temperature is found to be substantially higher 
than the currently accepted values, also obtained using 
QMC, but with a different data analysis; as expected, 
they are also substantially lower than estimates obtained 
within a Hartree-Fock/mean field approximation. 
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